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ABSTRACT 

In recent years it has become increasingly clear that Active Galactic Nuclei, and radio-galaxies in 
particular, have an impact on large scale structure and galaxy formation. In principle, radio-galaxies 
are energetic enough to halt the cooling of the virialized intracluster medium (ICM) in the inner re- 
gions of galaxy clusters, solving the cooling flow problem and explaining the high-mass truncation of 
the galaxy luminosity function. We explore this process through a series of high resolution, three di- 
mensional hydrodynamic simulations of jetted active galaxies that act in response to cooling-mediated 
accretion of an ICM atmosphere. We find that our models are incapable of producing a long term 
balance of heating and cooling; catastrophic cooling can be delayed by the jet action but inevitably 
takes hold. At the heart of the failure of these models is the formation of a low density channel 
through which the jet can freely flow, carrying its energy out of the cooling core. It is possible that 
this failure is due to an over-simplified treatment of the fast jet (which may underestimate the "dentist 
drill" effect). However, it seems likely that additional complexity (large-angle jet precession or ICM 
turbulence) or additional physics (magnetohydrodynamic effects and plasma transport processes) is 
required to produce a spatial distribution of jet heating that can prevent catastrophic cooling. This 
work also underscores the importance of including jet dynamics in any feedback model as opposed to 
the isotropically inflated bubble approach taken in some previous works. 

Subject headings: cooling flows — galaxies: active — galaxies: formation — galaxies: jets — hydro- 
dynamics — X-rays: galaxies: clusters 
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1. INTRODUCTION 

In recent years, it has become increasingly clear that 
Active Galactic Nuclei (AGN) may well have an impact 
on large scale structure and galaxy formation. The ac- 
cepted framework for galaxy formation states that bary- 
onic matter falls into developing dark matter halos and 
the resulting accretion shocks raising it to the approx- 
imately the virial temperature of the halo. The "cold" 
baryonic components of the galaxy then form via radia- 
tive cooling of this shocked gas. In the absence of any 
feedback processes (i.e., with gravitational collapse fol- 
lowed by radiative cooling alone), the galaxy mass func- 
tion would have to essentially follow the dark matter halo 
mass function — in essence, the baryons that are within 
the turn-around radius of the developing dark matter 
halo are trapped and fated to eventually form the bary- 
onic galaxy at the halo's center. However, as has been 
long known, this is clearly not the case. At both the 
highest and lowest masses, there is a deficit of galaxies 
compared with dark matter halos. These major discrep- 
ancies indicates that multiple, non-gravitational, feed- 
back mechanisms act on galaxy formation. 

At the low luminosity/mass end of the galaxy distri- 
bution, the deficit ca n be explained entirely by feedbac k 
from star formation l|Larsonl fl974t iDekel fc SiHd [l986L 
Once the first (massive) stars to form from the cooling 
baryons start to supernova, superwinds from the proto- 
galaxy can expel a significant fraction of the remaining 
baryons from the (shallow) gravitational potential well 
of the dark matter halo. The explanation for the deficit 
at the h igh mass end is not so clear, though. Inefficient 
cooling C Rees fc Ostrikerll977j) of baryons that have been 
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shocked to the high virial temperature (exceeding 10 7 K) 
of large dark matter halos can account for some of the 
deficiency but is insufficient. An additional mechanism 
is required that is substantially mo re efficient than sta r 
formation (e.g., see discussion in iBenson et alJ f2003). 
Many authors have suggested that radio-loud AGN pro- 
vide this additional feedback. 

If radio-loud AGN do indeed regulate the formation 
of the most massive galaxies, this is probably one of 
the few aspects of galaxy formation that can be stud- 
ied in detail in the local universe: i.e., in the cores 
of cooling galaxy clusters. Galaxy clusters possess the 
largest dark matter halos found in the Universe and con- 
tain virial baryons (the intracluster medium; ICM) at 
X-ray emitting temperatures. The ICM of galaxy clus- 
ters has been intensively studied by every X-ray observa- 
tory since its discovery in th e earliest days of X-ray as- 
tronomy ijFelten et al.l ll966). but imaging spectroscopy 
by Chandra and XMM-Newton have raised these studies 
to an unprecedented level. At its simplest, these obser- 
vations provide measurements of the density and tem- 
perature of the ICM, which allows a radiative cooling 
time to be computed. In the central regions of most 
relaxed clusters, the cooling time is often significantly 
shorter than a Hubble time (and often as low as few x I0 8 
years). With such short cooling times, there must be 
either a growing "sink" of cooled gas or a heat source 
acting on the cluster center. Observationally, there are 
very strong limits of the amount of cool gas present. 
In particular, there is usually insufficient star formation 
seen in the central cD galaxy to account for the cooling 
flow [although ther e is clearly star forma ti on oc curring; 
IQ'Dea et abl fffinks fc Mushotzkvl IpOfil) ]. The 

detailed discrepancy between the X-ray measured cooling 
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rate and the lac k of cooled ga s forms the classic "cooling 
flow problem" i|Fabiar]ll994fl . Dispersive spectroscopy 
by XMM-Newton deepened this mystery. While most 
clusters clearly show evidence for radiative cooling of the 
ICM from the virial temperature T V ir to about T v j r /3, the 
absence of lower-ionization iron lines sets tight limits on 
the emission measure of gas below T v j r /3 (which would b e 
expected to be a pro lific radiator; iPeterson et alJ l)2001|) : 
iTamura et all l|200im . 

Of course, there is an obvious link between the appar- 
ent inability of the ICM to cool below X-ray emitting 
temperatures and the deficit of massive galaxies. It is 
precisely the cD galaxies at the center of rich galaxy clus- 
ters which are "trying" to grow above the observed cutoff 
in the galaxy mass function via the radiative cooling of 
the core regions of the ICM. This connection strongly 
suggests that solutions to the cooling flow problem that 
rely on "hiding" the signatures of cooling (e.g., through 
strong metallicity inhomogeneities or mixing layers) are 
doomed to fail, thereby increasing the motivation for pro- 
viding a viable heat source for the cluster core. 

Radio-loud AGN are very energetic and are fre- 
quently found in the giant elli ptical (cD) galaxies i n 
the center of cooling clusters l|Burns fc Owenl Il977l) . 
which makes them ideal candidates for the heat source 
needed to offset cooling cluster. Even putting galaxy 
formation and cooling flow arguments aside, however, 
it is clear that radio loud AGN have an impact on 
their environm ent. First seen i n Ein st ein and RO SAT 
observations llFeigelson et alJ 119871 iBohrin ger et alJ 
Il993l H99l ICarilli et all 119941: IHeinz et alJ 11998). 
Chandra has studied nume r ous examples of ICM 
bubbles llFabian et al] l200(Tt IMcNamara et all 1200(1 
iBlanton et alJ 120011: lYouirg et al.ll2002D. ghost bubbles 
(IMcNamara et alJ 120011: IHeinz et all 120021 iChoi et al l 
120041 iFabian et all 12000). ripples llFabian et alJ 120031 
2005rJ~ shells l|Fabian et all 12000^ and filaments 
that are clearly associated with a central radio-loud 
AGN. There is a large body of numerical work which 
models various aspects of these interactions. Early 
hydrodynamic models of the putative jet- cocoon/shock 
struc ture in Cygnus-A were presented by ijClarke et al.1 
1997). More recent numerical investigations have 
explored the buoyant evolution of the cocoon af- 
ter i ts supersonic expansion phase ||Chunizov^^d 
20011 [Briiggen fc Kaiserl l200ll l2002t iRevnolds et al 
20021 IBasson fc Alexander! I2003t lOmma et al.l \2~ 



Omma fc Binneyl 12004 iDalla Vecchia et, alJ 20041 
Robinson et all 120041: iZanni et all 120051) . the ef f ect of 



lasma transport pr ocesses ijRnszkowski et akl 120041: 
Re^noJo^s_e^ li alJj20 i i 5H) and the action of magne tic fields 
(Rob inson et alJl2004t Uones fc De Youndl2005|) . 

Despite the reasons to look to AGN heating as the 
solution to the cooling flow problem (and hence the reg- 
ulator of high mass galaxy formation) , there is a problem 
with this picture. Direct observational signatures of ICM 
heating remain elusive in the observations despite the 
existence of some very deep Chandra and XMM-Newton 
observations of nearby vigorous ICM/ AGN interactions. 
Strong shocks are not seen in most systems even though 
they are seen in many simul ations that capture strong 
AGN heating (however, see lOmma et alJ |2004) for a 
discussion of the difficulty of seeing shocks in slow jet 
simulations). Furthermore, the recent sample of radio- 



galaxy blown ICM bubbles bv iBirzan et alJ l|2004|) sug- 
gests that the total work done by the AGN to inflate the 
observed ICM cavities may not be sufficient (by approx- 
imately an order of magnitude) to offset cooling in many 
clusters. This may point to the i mportance of heating 
by the dissipation of sound waves llFabian et al.ir2 005a) 
or the decay of globa l ICM mo des l|OimTraetnilTl20C3l 
IRevnolds etaI1l2005a|) 

In this paper we perform high-resolution, three di- 
mensional ideal hydrodynamic simulations of AGN feed- 
back in a relaxed cooling cluster. Unlike many of 
the previou s 3-d investigations [ with t he notable ex- 
ceptio n s of IBass on fc A lexander! l)2003fl : lOmma et alJ 
(120041k lO mma fc BirmevT l)2004|) ] . we simulate the injec- 
tion of a supersonic jet into the ICM atmosphere rather 
than starting with an initial condition of a pre-inflated, 
static radio cocoon. We consider this an important issue 
— only through a direct modeling of the jet can we hope 
to be able to capture the jet-induced ICM shock heating 
as well as the complex internal dynamics of the cocoon. 
While the effect of the jet on the ICM is treated from 
first principles (through the evolution of the equations of 
hydrodynamics), the enormous range of scales between 
the ICM core and the gravitational sphere of influence of 
the black hole prevents us from treating the formation 
of the jet or the AGN fueling properly. In this paper, 
we introduce several different feedback prescriptions by 
which the AGN jet power is related to the ICM proper- 
ties in the innermost region. Our goal is to construct a 
model system in which AGN heating and ICM cooling 
are, in the long term, in balance. We find that our ideal 
hydrodynamic models fail to achieve this balance. As- 
suming that the AGN feedback picture is indeed correct, 
we conclude that our models must be failing to capture 
some crucial aspect of the system. 

In Section |2 the basic setup for our simulations is de- 
scribed, including a description of the model ICM atmo- 
sphere and our hydrodynamic code. Section 3 presents 
our results for the different feedback scenarios consid- 
ered. We discuss our results along with possible ways to 
rescue the AGN feedback hypothesis in Section^] Finally 
our modifications to the ZEUS-MP code are explained in 
Appendix lAl 

2. BASIC SETUP 

We aim to model the ICM of a cluster and its inter- 
action with a central radio galaxy. Our basic picture is 
a spherical cluster consisting initially of stationary dark 
matter and gas. The gas is initially in hydrostatic equi- 
librium, but is cooling through optically thin thermal 
emission. As the gas cools, some of it will flow across the 
inner radial boundary of the simulation and is no longer 
simulated. The amount of gas to cross this boundary is 
used as our primary diagnostic of the cooling flow. 

Initially, the cluster is spherically symmetric and 
isothermal. The gas is setup with a /3-model profile, 
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This is a simple analytic fit to cluster mass profiles. The 
gravitational potential, 
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is the set so the initial gas configuration is in hydrostatic 
equilibrium. The potential is assumed to be generated by 
a static distribution of dark matter which remains fixed 
throughout the simulation. The self gravity of the gas 
is ignored (we note that in a standard cosmology, the 

;as mass is only aro und 14% of the total cluster mass 

Ostriker et"al]l2005|l ). 

Spherical polar (r, 8, 4>) coordinates are used. All sim- 
ulations were run on a 200 x 200 x 100 cell grid, with 
enhanced resolution near the center of the cluster and 
near the jet axis. This corresponds to a physical grid 
with a inner radius at r — 10 kpc and an outer radius 
at r = 1000 kpc. The angular coordinate was only al- 
lowed to vary from 9 = to 9 = §. This effectively only 
covers half the cluster and therefore only allows for one 
jet. Reflecting boundaries were use to mimic the effect of 
the missing half of the cluster. This aids in maintaining 
a high-resolution while keeping a reasonable number of 
grid cells. 

In ideal (adiabatic) hydrodynamics, the equations of 
hydrodynamics can be written in dimensionless form 
and, hence, the results of a single simulation may be 
scaled to a whole family of prob lems with different mas s, 
size and time scales (e.g., see iRevnolds et aD (J2002)). 
Since we must add radiative cooling to our simulations 
(see Section 13.10 , we are forced into a set of physical 
units. The values for the rich cluster of [Reynolds et al.l 
( 2002) were used. This gives us a core radius of r D = 
100 kpc, a code length unit of r = 50 kpc, time units 
of 50Myrs, a sound speed c s = 1000 km s -1 . The re- 
sulting total radiative luminosity of the model ICM is 
1.22 x lQ 45 ergs _1 . 

2.1. Code 

The simulations were performed using the ZEUS-MP 
code. A modified version of the original National Cen- 
ter for Supercomputing Applications (NCSA) release was 
used. We have made our modified version publicly avail- 
able 1 . A further discussion of this version of ZEUS-MP 
detailing our modifications is given in Appendix 1X1 

ZEUS-MP is a parall el version of the ZEUS m agneto- 
hydrodynamic code l|Stone fc Normanlll992alH) . ZEUS 
is a fixed-grid, time-explicit Eulerian code which uses an 
artificial viscosity to handle shocks. When operated us- 
ing van Leer advection, as in this work, it is formally of 
second order spatial accuracy. The work reported here 
used this code in a pure hydrodynamic mode. All simu- 
lations were done on the University of Maryland, Astron- 
omy Department's GNU/Linux Beowulf cluster with run 
times of the order of one to two months for four processor 
simulations. 

ZEUS solves the standard hydrodynamic equations, 
which, including radiative cooling fSection l3.ll) are 
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Name Feedback 


Efficiency Delay 


Rotation 


Run A dud i 


NA 


NA 


NA 


Run B strong jet 


NA 


NA 


NA 


Run C weak jet 


NA 


NA 


NA 


Run D simple feedback 


0.0001 


NA 


NA 


Run E simple feedback 


0.00001 


NA 


NA 


Run F delayed feedback 


0.0001 


short 


NA 


Run G delayed feedback 


0.00001 


short 


NA 


Run H delayed feedback 


0.0001 


long 


NA 


Run I delayed feedback 


0.00001 


long 


NA 


Run J delayed feedback 


0.01 


long 


NA 


Run K delayed feedback 


0.1 


long 


NA 


Run L delayed feedback 


0.00001 


long 


solid body rotation 


Run M delayed feedback 


0.00001 


long 


Fraction of grav. rotation 



where 
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dt 
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Stability is determined by the usual Courant-Friedrichs- 
Lewy (CFL) condition. As the jet is typically very su- 
personic, its speed provides the limiting timestep during 
most of the simulations. 

3. SPECIFIC MODELS 

A set of thirteen 3-d simulations were performed. Ta- 
ble n lists the simulations along with some relevant pa- 
rameters. In the following sections we will describe each 
simulation along with presenting the results of each. 

An important diagnostic used to compare the differ- 
ent models was the mass accretion rate across the inner 
boundary of the simulated grid. This was measured in 
M Q year _1 . To calculate the mass flow, the amount of 
mass in each cell on the inner boundary with a negative 
(inward) velocity in the radial direction was summed. 
This is probably not completely accurate at the highest 
densities and mass accretion rates as the inflow bound- 
aries are not perfectly efficient and sound wave may be re- 
flected off the boundary. However, this should not change 
the results since these inaccuracies do not occur until the 
mass flow and density reach physically unrealistic values. 

3.1. Radiative cooling 

The driving force behind the cooling flow is the thermal 
bremsstrahlung and line radiation which removes ther- 
mal energy from the ICM core. This is modeled with the 
optically thin cooling law 

A = [Ci(fe B r) Q +C 2 {k B Tf + C 3 ]0.704 {p/m p f 

xl0~ 22 ergs cm~ 3 s _1 , (7) 

which is the same law used bv iRuszkowski fc Beeelmanl 
(2002) with the coefficients d = 8.6 x 10" 3 , C 2 = 
5.8 x 10~ 2 , C 3 = 6.4 x 1Q- 2 , a = -1.7, and = 0.5 
({Sutherland k, Dopital ll993'). This cooling term is pro- 
portional to p 2 which produces the eventual catastrophic 
cooling as the center of the cluster becomes denser as it 
cools. Below a minimum temperature (0.1 keV), cooling 
was manually truncated. Even if we allowed material to 
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cool below this limit, our spatial resolution would be in- 
sufficient to follow the resulting structures and additional 
physics not captured by eqn.[7| would become applicable. 

The cooling represents an extra term in the energy 
equation (eqn. EJ) and is implemented as an explicit 
source term. As with all physical processes, there is a 
maximum allowed time step for numerical stability asso- 
ciated with the cooling, but in our case it is always above 
the normal hydrodynamic CFL condition and need never 
be implemented. 

Due to the density dependence of the radiative cooling, 
it is only relevant in the inner regions of the cluster. This 
is important for three reasons. First, as the core cools, 
the inner regions loses pressure support. This causes the 
inward flow of material, some part of which could fuel 
the accretion on the central compact object. Second, 
as the inner regions become denser due to the sagging, 
the cooling increases. This causes the eventual runaway 
cooling whose absence represents part of the cooling flow 
problem. Finally, since the outer region does not cool 
appreciably with a Hubble time or so, it represents are 
large reservoir of hot material which could potentially be 
exploited to balance cooling (although this does not seem 
possible in purely hydrodynamic models). 

To establish a control, our first simulation (Run A) 
followed the pure radiative collapse of our model ICM 
atmosphere (i.e., a spherically symmetric homogeneous 
cooling flow). The evolution of this system was very sim- 
ple. The cluster cools, primarily in the inner regions. As 
it cools, the inner regions become denser. As the gas be- 
comes denser, it cools quicker, and the process runs away 
(until numerical issues with the very cool dense gas force 
us to terminate the simulation) . The mass accretion rate 
for pure cooling can be seen in Figure ^ With no mech- 
anism to stop or slow the cooling, the mass accretion 
show a featureless, approximately exponential increase. 
This unbounded cooling and mass accretion is our first 
example of a catastrophic cooling. 

To diagnose the time taken for the system to undergo a 
cooling catastrophe, we measure the time taken for each 
of our simulations to exceed a mass accretion rate of 
SOOOMoyr^ 1 (an arbitrary "large" mass accretion rate). 
For the pure cooling flow run (Run A), this threshold 
mass accretion rate is crossed at 244.1 Myrs. Mate- 
rial first fell below our imposed lower temperature limit 
at 248.9 Myrs. The radial dependence of temperature 
for Run A is shown in Figure Initially the cluster is 
isothermal. The inner regions cool first, while the outer 
regions barely change in temperature. As catastrophic 
cooling occurs, the temperature gradient in the inner re- 
gion becomes progressively steeper. 

3.2. Single AGN Outburst Models 

The onset of the cooling catastrophe in the homoge- 
neous cooling flow is well known and unsurprising. Our 
hypothesis and the focus of this work is the idea that jet 
activity by a central AGN can heat the ICM core and pre- 
vent this catastrophic cooling. As an initial exploration 
of jet heating models, we follow the evolution and effects 
of a period of jet activity in which the jet has a fixed and 
constant power for a pre-defined duration. In Run B, the 
jet is active for 50 Myrs after which it is completely shut 
off and the resulting ICM allowed to evolve passively. 
These cases are essentially 3-dimensional generalizations 



of the axisymmetric simulations of lRevnolds et all l)2002fl 
and share ma ny characteristics with (but are higher res- 
olution than) iBasson fc Alexa nder! 12 003) . The axisym- 
metric simulations of lZanni et alJ l)2005[) provide another 
example of this type of single burst jet heating. This can 
be viewed as a zeroth order approximation to AGN feed- 
back. 

In detail, the action of the AGN is modeled as a jet of 
low density material injected from the inner radial edge 
of the simulation grid using an inflow boundary condi- 
tion. The jet has a density of 1/100 of the initial inner 
ICM density and (at the inner boundary) is in pressure 
equilibrium with the initial inner ICM. Therefore, the in- 
ternal sound speed of the injected jet material is 10 times 
that in the initial ICM atmosphere. The jet has an open- 
ing angle of 15°. In the case of a single jet outburst of 
Run B, a Mach number of around 10.5 compared to the 
background material was used resulting in a total kinetic 
luminosity of 9.8 x 10 45 erg s _1 . We also modeled a rather 
weaker jet with kinetic luminosity of 8.9 x 10 45 ergs -1 in 
Run C (which was otherwise identical to Run B). 

Plots of entropy are show n at various times for Run B 
in Figure As discussed in lRevnolds et al.1 l|2002j) . plots 
of entropy are a good tool for differentiating shocked 
jet material from background material. Initially, the 
jet carves a path through the ICM, and terminates in 
a shock. It can be seen that the jet channel is sur- 
rounded by a backflow of shocked (high entropy) jet ma- 
terial. This shocked material is over-pressurized com- 
pared to the background, and expands laterally into a 
cocoon structure. The cocoon is separated from the back- 
ground by a contact discontinuity. Rayleigh- Taylor (RT) 
and Kelvin-Helmholtz (KH) instabilities work to shred 
this cocoon and mix in background material. After the 
jet has been turned off, the cocoon is left to evolve pas- 
sively. Buoyancy forces cause it to rise, leaving the clus- 
ter core. Once the cocoon completely detaches itself from 
the core, it takes on the appearance of a rising bubble, 
which spreads out and fades somewhat into the cluster 
background (much like observed ghost bubbles) until it 
leaves the computational grid. 

Figure 0] shows mass accretion for Run B and Run 
C. Although these differs slightly in the timing, they 
qualitatively show the same behavior. Initially, there is 
moderate mass flow (a few hundred M Q yr _1 ). Ironically, 
the jet-activity is actually responsible for this initial en- 
hanced period of inflow; some matter in the innermost 
region of the ICM core becomes caught in the backflow 
that results from the onset of the of the jet activity and is 
swept across the inner boundary. After the initial spike, 
the mass accretion rate drops off to very low values while 
the jet is on and remains very low for over 100 Myrs after 
the jet stops due to the shock heating of the ICM core. 
The cooling flow then starts to re-establish and the mass 
flow begins to increase (eventually) catastrophically as it 
does in the pure cooling case. The time for catastrophic 
cooling has been delayed somewhat from the pure cooling 
case. For Run B, M does not reach 5000M Q yr _1 until 
306 Myrs (more than 60 Myrs later than the pure cooling 
model). Material falls below the minimum cooling tem- 
perature at 305 Myrs (i.e., at essentially the same time 
that the mass accretion rate crosses our "catastrophe" 
threshold). The weaker jet, Run C, reaches the catas- 
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Fig. 2. — Radial temperature dependence for pure cooling (Run A). 



trophic point at 285 Myrs, which is only slightly delayed 
from the pure cooling case. The temperature falls below 
our floor at 274 Myrs (also slightly later than the pure 
cooling case). 

3.3. Feedback Models 

The radio-loud AGN feedback hypothesis would ar- 
gue that the central AGN acts such as to maintain a 
long term balance of jet heating and radiative cooling of 
the ICM. Furthermore, this must occur in systems with 
widely varying radiative luminosities. Clearly, this mech- 
anism requires that the AGN power be somehow regu- 
lated by conditions within the ICM so that it does not 
underheat or overheat the ICM core [see lBinnev fc Taborl 



l|1995(l for early work on AGN feedback on cooling flows] . 
In this work, we explore a set of models in which the 
kinetic luminosity of the jets is connected to the mass 
accretion rate in the cooling flow. Unlike the other ef- 
fects (cooling, jet propagation, etc.), feedback is the only 
one that does not directly represent fundamental physics 
of the cluster gas. Instead, it tells us something about 
fueling in the central engine of the AGN. This provides 
us with the most freedom in how to implement it. Due 
to the large range of scales involved, it is not possible 
to directly model feedback. This would require resolv- 
ing gas for the outer edges of the cluster (nearly a mega- 
parsec) all the way down to the inner accretion disc (par- 
sec scales). A direct assault on this problem would not 
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be practical at the current time even for an Adaptive 
Mesh (AMR) code. Instead, we must use some simpli- 
fied model for feedback. 

In our first set of feedback simulations (Runs D and 
E), the jet is injected with a kinetic energy proportional 
to the instantaneous mass accretion rate across the in- 
ner radial boundary of our simulation domain. This is 
achieved by modulating the injection velocity to be 
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2r]Mc 2 



(8) 



where M is the mass flow rate across the inner radial 
boundary of the simulation, c is the speed of light, p is the 



density of the injected jet material, A is the area of the 
jet "nozzle" on the inner boundary, and rj is the efficiency 
with which the rest mass energy of the ICM cooling flow 
is converted into jet power. The appropriate choice for rj 
is far from clear. In the extreme case that all of the mat- 
ter that flows across the inner boundary was to accrete 
onto the central supermassive black hole, r\ would be the 
jet-production efficiency of the actual black hole accre- 
tion disk itself. For an efficient disk, we would have rj ~ 
0.1, although hi gher efficiencies are possi ble if the black 
hole is spinning l|Novikov fc Thorneill97^ or if magnetic 
coupling wi thin the radius of marginal stability become s 
important l|Gammidll999t lArmitaee fc Reynolds! 12003). 
However, as has been extensively discussed, supermas- 
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sive black holes in the centers of giant elliptical galax- 
ies accreting at a modest rate from the hot interstellar 
medium may well have a significantly smaller efficiency 
(j) ~ 10 -3 ). Much of the remaining energy is advected 
across the event horizon or driving a slow, uncollimated 
wind from the disk. Furthermore, it is plausible that 
only a small fraction of the mass that flows across the 
inner boundary of our simulation domain (located at 
10 kpc) actually enters the sphere of influence of the mas- 
sive black hole, with some of the remaining going to fuel 
low-level star formation. 

Motivated by this discussion, we explored two cases 
of simple feedback with efficiencies of r] = 10 -4 (Run 
D) and r\ = 10~ 5 (Run E). The resulting mass accre- 
tion for Run E (lower efficiency) is given is Figure [SJ 
The higher efficiency (Run D) shows the same pattern. 
Initially, the mass accretion rate cycles up and down in 
response to the jet power (partly due to the backflow 
noted above), as we would expect if feedback works to 
stop cooling. However, by 150 Myrs, this breaks down. 
From that point on, the mass accretion rate increases, 
and the increasingly powerful jet does nothing to stop it. 
The time for catastrophic cooling and for material to fall 
below the minimum cooling temperature is not signifi- 
cantly affected by this type of feedback. Both runs have 
very similar catastrophic cooling times with 246.7 Myrs 
for Run D and 245.43 for Run E. This is essentially the 
same as the time for the pure cooling model. 

The reason for this failure of feedback to prevent the 
cooling catastrophe can be seen in Figure H3 After the 
jet has been active for a while (even at very low power) it 
clears a low-density channel in the ICM. As cooling and 
ICM accretion proceeds along equatorial latitudes (i.e., 
in the plane perpendicular to the jet axis), the increas- 
ingly powerful jet flows freely down the pre-cleared chan- 
nel. This prevents the jet from depositing energy near 
the cluster core. Instead the kinetic energy is carried to 
the head of the cocoon reasonably unimpeded. With no 
energy deposition near the core, cooling proceeds on the 
catastrophic course in the equatorial regions almost as if 
there were no jet. 

3.4. Delayed Feedback Models 

The failure of the simple feedback model is, at least in 
part, related to the formation of t he low density chan- 
nel along the jet axis. We note that lOmma et all ( 2004) 
also discuss the importance of such channels within the 
context of their slow-jet simulations and argue that the 
effective heating required the channels to fill in between 
powerful outbursts. Motivated by this, and in an at- 
tempt to model more realistic feedback, we performed a 
set of simulations in which a time delay was introduced 
between the mass accretion rate and the response of the 
jet. Clearly the immediate feedback of Run D and E is 
not physically accurate. Some time delay must be added 
to account for the material travel time from the inner 
edge of the grid to the accretion disk and onto the black 
hole. We note that the time for the relativistic jet to 
reach from the black hole and to enter the computational 
domain is negligible. 

In Runs F and G, a delay of 10 Myrs was introduced. 
This is the sound crossing time of the cluster center (the 
area from the inner edge to r = 0) and represents the 
minimum physically reasonable delay. Runs H, I, J, and 



K all have a delay time of 100 Myrs. This is approxi- 
mately the dynamical time for the central galaxy and is 
a reasonable "best-guess" at a physically plausible time 
delay. Interestingly, this is also the approximate cool- 
ing time observed by Chandra in the centralmost regions 
of cluster cooling cores. The efficiencies of these runs is 
listed in Table □ 

An example of the low efficiency delayed feedback with 
a short delay can be seen in Figure [7| There is not 
much difference between this and the immediate feed- 
back. There is also no essential difference in the time to 
reach catastrophic cooling. The only minor effect we see 
here is that the lower efficiency run (G) has a slightly 
lower rate of increase for the accretion, probably due to 
less efficient channel formation. 

For Runs H through K, the physically motivated de- 
lay of 100 Myrs was used. Runs H and I and essentially 
copies of Runs F and G with a longer delay between the 
mass flow and the jet. As Figure |H1 shows, this have a 
more significant effect on the mass flow than the shorter 
delay had. For the first 100 Myrs, there was no jet, so 
the run proceeded as a pure cooling run. The onset of 
(weak) jet activity caused a sharp dip in the accretion, 
followed by a spike (up to 600 M Q year _1 ) lasting for 
about 50 Myrs (or half of the dynamic time). However, 
after several smaller smaller spikes and dips, runaway 
cooling proceeds in the equatorial plane. In these cases, 
the cooling catastrophe is delayed as compared with the 
pure cooling model, but only by 40 Myrs. The higher effi- 
ciency version (Run H) behaved similarly only it reached 
the cooling catastrophe slightly faster (once again due to 
increased channel formation). 

Given the difficulties with the feedback models found 
above, we also ran two simulations with very efficient 
jet production; Run J and K possessed rj = 0.01 and 
0.1 respectively. Run K was performed in the spirit of 
an extreme case and we note that it is not reasonable 
that accretion onto the central galaxy and then onto the 
central AGN are both perfect, with the only loss coming 
from relativity and the central black hole. Even at an 
efficiency of 0.01 it is hard to come up with a plausible 
scenario for such nearly perfect accretion. 

Even the extreme runs are not able to completely stop 
cooling. For Run J, the temperature falls below T m - m at 
245 Myrs. This is even sooner than in the pure cooling 
case. Run K does not have this until 267 Myrs which is 
between the times for the other delayed runs and the sin- 
gle jet runs. The extreme efficiency feedback does, how- 
ever, significantly delay catastrophic cooling. Run J does 
not experience catastrophic cooling until 345 Myrs. The 
mass accretion behavior is shown if FigureEl Although it 
has taken longer for the cooling to reach our catastrophic 
level, the final catastrophe is extremely sharp, with the 
mass accretion rate exponentiating on timescales of only 
1 Myr or so. This is much more dramatic than any of 
the less efficient runs. Run K holds of catastrophic cool- 
ing for 50 Myrs beyond Run J (398 Myrs), but it too 
becomes catastrophic rapidly at that point. 

3.5. Feedback with Rotation 

The AGN feedback and jet heating models described 
above clearly fail to adequately heat the ICM core in 
a manner that stably offsets radiative cooling. Our fi- 
nal two simulations (Runs L and M) explore the effect 
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Fig. 5. — Mass accretion rate for feedback model (Run E). 
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Fig. 6. — Temperature Map (top), Pressure Map (middle), Entropy Map (bottom) for immediate feedback jet. The temperature map 
only shows the highest temperatures to pick out primarily jet material. The thin, low density channel can be easily seen in temperature 
and pressure (and to a lesser extent entropy). Only the inner 254 kpc of the simulation is shown. 



of an ICM atmosphere that is not initially static. The 
ICM of real clusters is never, of course, a perfect hydro- 
static atmosphere. Clusters have complex merger and 
formation histories that leave an imprint on the dynam- 



ics of the ICM {B urns &: Owenlll977|) . In several clusters 
observed with Chandra, the ghost cavities are not svm- 
metric about the c luster center (e.g., A 4059, lHeinz et alJ 
( 2002); Perseus- A. iFabian et alJl)2003|) ). suggesting that 
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Fig. 8. — Mass accretion rate for low efficiency long delayed feedback (Run I). 



they are buoyantly rising in an ICM atmosphere which 
itself has velocity structure. While a full exploration of 
this class of models is beyond the scope of this paper, our 
final two simulations model the case of a rotating ICM 
atmosphere. Due to the assumed reflection symmetry in 
the plane perpendicular to the jet axis, we are forced to 
consider only rotations with an axis coincident with the 
jet axis. 

In Run L, a delayed feedback model similar to Run 
I is examined in which the initial ICM atmosphere is 
assumed to undergo solid body rotation. The rotation 
speed at the outer edge of the grid is set equal to the 
sound speed of the cluster. This leave the rotation in 
the core at a very low value. The results from this were 



nearly indistinguishable (in cooling, not in detailed phys- 
ical structure) from Run I. This is because the rotation in 
the core was not high enough to mix the ICM effectively 
or stop accretion in any other way. Without having gas 
moving supersonically in the outer parts of the cluster 
(which is clearly unphysical) , there is no way to get solid 
body rotation to help offset cooling. 

In Run M a rotation law was chosen such that the 
centrifugal force associated with the rotation is a fixed 
fraction (10%) of the gravitational force at that loca- 
tion. This results in a rotation law in which the angular 
velocity increases with decreasing radius until one gets 
well within the ICM core, at which point the angular 
velocity smoothly goes to zero at the center. Figure ITUl 
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Fig. 9. — Mass accretion rate for very high efficient feedback (Run J). 



from Run M shows the first accretion curve that differs 
significantly from the previous patterns. Initially, the ac- 
cretion rate goes up and down in response to the jet. At 
around 240 Myrs, gas falls below the cooling limit (sim- 
ilar time to the pure cooling and most feedback runs). 
At 285 Myrs, the accretion passes the 5000 M Q year _1 
limit and continues to climb. Unlike previous setups, it 
does not grow without bound (or flatten out at some un- 
realistically high value). Instead at 16,000 M©yr _1 , the 
curve reverses and over the course of 50 Myrs drops back 
down to very low levels (a few hundreds of solar masses 
per year. It then continues at that rate with only small 
fluctuations for the duration of the simulation. 

In this case, the AGN behavior is actually rather inci- 
dental. At early times, the slowly rotating inner core of 
the cluster cools and accretes. Higher angular momen- 
tum material flows inwards and, eventually, the ICM core 
becomes rotationally supported. Given the assumptions 
of our model, this material will conserve its angular mo- 
mentum and, irrespective of radiative cooling, will not 
accrete. The main role of the AGN is to prevent accre- 
tion of ICM from the high latitude regions within the cen- 
trifugal barrier. Indeed, the cooling catastrophe has not 
in fact been averted. Rapid cooling into a rotationally 
supported disk can clearly be seen in the late stages of 
this simulation fFigure lTT|) . Our mass accretion diagnos- 
tic (based on mass flow across the inner radial boundary 
of the simulation) fails to detect this particular manifes- 
tation of the cooling catastrophe. 

4. DISCUSSION 

Our work with simple hydrodynamic models show 
that it is more difficult than previously assumed to halt 
a cooling flow through AGN activity. The ultimate 
h ope (as in, for example, t he on e dimensional models 
of Ruszkowski & Bcgclman (2002)) is that a proper cou- 
pling of the AGN kinetic luminosity to the cluster gas 
will allow a long-term balance to be established where 
the AGN heating balances the ICM radiative cooling. In 



our three dimensional simulations, we do not see such 
balance. Instead, the AGN jet seems to effectively clear 
a low density channel through the ICM core, after which 
time the kinetic energy of the jet is deposited well outside 
of the cooling core. Thus, simulations with delayed feed- 
back do not produce many multiple outbursts. Instead 
they produce one or two long outbursts with varying jet 
power before the ICM core catastrophically cools. 

Our modeled AGN feedback is not totally ineffectual 
at heating the ICM. Both the time for catastrophic cool- 
ing and the time for gas to fall below our cooling limit 
are increased by the AGN in most cases. The single jet 
cases are rather effective (delaying the catastrophe for 
as long the most effective feedback simulations). How- 
ever, the instantaneous feedback and the feedback with 
a time delay of the core sound crossing time have ba- 
sically no effect of the catastrophe time (no more that 
a few Myrs). The more realistic delays are capable of 
delaying the cooling catastrophe by several tens of Myrs 
(more in the unrealistically efficient cases). The only ex- 
ception to this behavior was the rotating cluster modeled 
in Run M. As was shown in Figure there is a brief 
period of very high accretion, following by an extended 
period of lower accretion. However, as discussed above, 
catastrophic cooling is still occurring into a rotationally 
supported disk. The rotational support prevents this gas 
from falling into the central region, but it is still a large 
body of cool and cooling gas around the cluster center 
which contradicts observations. 

Due to the computationally expensive nature of our 
simulations and the multiple parameters present, it was 
not possible to test every permutation of the model pa- 
rameters. It is possible that there is some preferred area 
of the parameter space where the AGN heating balances 
the radiative cooling without producing unrealistic jets. 
The potential existence of such a privileged area of the 
parameter space does not change our conclusion as it 
would bring up a serious fine tuning problem. Cooling 
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Fig. 10. — Mass accretion rate for rotating cluster with feedback (Run M). 




Fig. 11. — Dense cold disk formed in the inner regions of Run M. 
is shown. 



Density at 880 Myrs shown. Only the inner 254 kpc of the simulation 



and AGN occur in a wide range of clusters of varying 
masses and temperature. For a regulatory mechanism to 
work, it must be general enough to work in the different 
clusters. We believe that our simulations sample physi- 
cally plausible parts of parameter space, thus the failure 
to model successful AGN feedback is an important result. 

As discussed in the Introduction, however, there is sig- 
nificant observational support for the notion that radio- 
galaxy heating can provide a long-term balance against 
radiative cooling. Thus, it seems clear that our models 
are not capturing some ingredients which are important 
in real systems. In the rest of this section, we explore 
some of these missing ingredients. 

We must acknowledge that while we do model the 
jet, we may well not capture all of the relevant jet dy- 



namics. Real AGN jets have internal dynamics that 
can only be captu red by considering special-relativistic 
I Marti et ail 19971) and magnetohydrodyna mics effects 
I Omma et al 1120041 lOmma fc Binnevl I2004J) . Further- 



more, our simulation probably has insufficient resolution 
to capture the full internal dynamics of a hypersonic jet 
even within the context of ideal non-relativistic hydro- 
dynamic. These effects may allow the jet momentum 
to be spread over a larger working surface (i.e., "dentist 
drill effect"), thereby helping alleviate the problems as- 
sociated with the jet drilling out of the ICM core. The 
fact that our simulations produce structures (cocoons 
and ICM cavities) with aspect ratios comparable to those 
observed, however, suggests that this is not a major de- 
ficiency of our simulations. 
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Before leaving issues associ ated with the prop erties 
of the jet itself, we not e that lOmma et alJ l)2004fl and 
lOmma fc Binnevl |2004) simulate the effects of slow (~ 
10, 000 kms -1 ), heavier jets on a cooling ICM core. In 
agreement with our conclusions, they also identify the 
importance of channel formation and the need for the 
channels to collapse between AGN outbursts. A new phe- 
nomena that appears to emerge from these simulations 
is the excitation of large-scale <?-modes in the ICM which 
persist (and slowly dissipate their energy in to the ICM) 
long after the AGN activity has shut-down. Tentative 
observational evidence for g-modes may exist in the clus- 
ter sur rounding the radio-galaxy 3C 401 l)Revnolds et alJ 
I2005al) . 

Any aspect of the system which alleviates the problem 
of channel formation will be potentially important in en- 
suring the jet energy is thcrmalizcd in the ICM core. The 
initial cluster in our simulations start with a smooth, hy- 
drostatic atmosphere. Real clusters have complex merger 
histories and turbulent ambient motions. This turbu- 
lence could help to distribute the energy from the jet 
more uniformly, or to oppose the formation of channels. 
Since the jets are powerful and higher velocity than even 
mildly supersonic ambient motions for most of their life- 
time, it is uncertain how effective turbulence would be 
in preventing channels and distributing energy. More 
importantly, real AGN jets may precess over significant 
an gles (for a recent cla im of jet precession in Perseus-A, 
see iDunn et al.l l|2006)). ICM channels may become ir- 
relevant and the kinetic energy thermalized in the ICM 
core if the jets can indeed precess over large angles on 
a timescale comparable to or less than the inner ICM 
cooling time. 

Finally, our results might be pointing to the impor- 
tance of physics beyond ideal hydrodynamics. A very 
real possibility is that plasma transport processes may 
be crucial. Thermal conduction provides a way to move 
thermal energy across a temperature gradient in the ICM 
and may be essential in helping the system avoid the 
cooling catastrophe. There is a large heat reservoir in 
the outer parts of the cluster which could be exploited 
to warm the core. Furthermore, thermal conduction and 
viscosity could be important mechanisms for thermaliz- 
ing the sound waves and weak shocks driven into the ICM 
by the radio-galaxy activity. All of these processes may 
delay the ICM cooling sufficiently to allow the jet-blown 
channels to refill, thus allowing a long term balance to be 
established. We note that there is growing observational 
evidence for the actio n of thermal conduction in real clus- 
ters (e.g., Perseus-A: iFabian et alJ l)2005blafl ). Magnetic 
fields are also certainly present in the ICM and, at some 
level, an MHD (or even a kinetic theory) treatment may 
be necessary. This will have an important effect on the 
stability of the buoyantly rising bubbles and the mix- 
ing of the high entropy radio plasma with the thermal 
ICM. But magnetic fields may have more profound con- 
sequences. Thermal conduction and plasma viscosity will 
both become very anisotropic due to the suppression of 
transport processes perpendicular to magnetic field lines, 
leading to qualitatively new flu i d inst abilities [the mag- 
netothermal instability Balbus ( 2000) and the fire-hose 
instability] . Cosmic ray pressure may also be relevant to 



driving ICM convection l)Chandranl f2004'). In a cooling 
cluster core with an embedded AGN, one can speculate 
that these instabilities drive ICM turbulence which may, 
itself, be a crucial ingredient in the dynamics of the sys- 
tem. Of course, this would be above and beyond any 
turbulence or bulk motions due to the merging history 
of the cluster. 

An important consequence of this work is that it firmly 
underscores the inadequacy of AGN feedback simulations 
that model AGN feedback through isotropically-inflated 
bubbles at prescribed locations in the ICM. Models that 
involve pre-inflated bubbles, or injecting energy isotrop- 
ically at given locations in the ICM, will fail to capture 
precisely that aspect of the physics that turn out to be 
crucial to the failure of our models — the development of 
a low density channel along which the AGN outflow can 
travel unimpeded thereby carrying its kinetic energy out 
of the cooling core. For this reason, it is vital that jet 
dynamics be included if a model is to properly address 
AGN-halted cooling in a cluster core. 

5. CONCLUSION 

We have performed a set of high-resolution three di- 
mensional simulations of jetted-AGN embedded in the 
cooling ICM cores of galaxy clusters in which the AGN 
power reacts in response to the cooling of the cluster 
gas. We make the first attempt to model regulated AGN 
feedback in the hope of achieving a long-term balance 
between radiative cooling and jet heating. However, we 
fail to construct a model in which AGN heating comes 
into long-term balance with the radiative cooling. The 
early time jet activity is extremely effective at clearing 
a channel through the ICM core. The jet at later times 
flows freely down this channel and hence deposits its ki- 
netic energy well outside the region which is starting to 
undergo catastrophic cooling. This essential behavior is 
robust to changing the efficiency or time delay character- 
izing the AGN feedback. Hence, we argue that some vital 
ingredient of AGN feedback is missing from our model; 
this may include more realistic jets, ambient ICM mo- 
tions/turbulence, jet-precession, or additional ICM mi- 
crophysics.Our results also highlight the absolute neces- 
sity to follow the jet dynamics in any attempt to address 
radio-galaxy heating of the ICM. 
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APPENDIX 
A. ZEUS-MP 

All simulations in this work were done using a modified version of the ZEUS-MP code. Our version is based 
on the initial NCSA released version (1.0b). Our modifications, documentation, and several supporting scripts have 
been made publicly available at http://www.astro.umd.edu/~vernaleo/zeusmp.html under the same terms as previous 
ZEUS releases. 

Here, we describe our changes and modifications to ZEUS-MP. First, the code had to be ported to compile with a 
current FORTRAN 77 compiler. As different compilers implement the standard (and the extensions to it) differently, 
we had to choose certain compilers as our targets. The Intel Compiler was used for the speed of the executables it 
produces for x86 compatible machines. For the sake of portability, we also maintain compatibility with the GNU 
compilers. None of these changes modify the behavior of the code. Primarily this involved removing multiply defined 
variables and cleaning up the namelist routines and the namelists themselves. Also, all filenames that differ only by 
case (common in the FORTRAN 77 build process) were changed to allow building on non-case sensitive filesystems. 

To allow for long runs, the restart routines (which did not work for parallel simulations) were completely replaced. 
The new restart routines work for parallel simulations and write to alternating files to save disk space. A wrapper 
script written in Perl is provided to correctly pick the most recent restart dump (if present) and start ZEUS-MP 
using that dump. The wrapper script also performs some basics checks of the integrity of the restart dumps before 
using them as incomplete or missing restart dumps have proven far more likely in parallel simulations than in single 
processor simulations. 

Despite the fact that ZEUS is written nearly entirely in FORTRAN 77, the C preprocessor is used heavily to allow for 
conditional inclusion of code. Some portions of this preprocessor code needed to be fixed; many preprocessor directives 
were not properly matched up or where improperly nested resulting in preprocessor flags that only worked for one of 
the two possible values. Also, some preprocessor directives did not enclose all of the code relevant to a given option. 
This made it impossible to completely turn certain options off in the code. The post-processor was modified to be 
useful for large number of output files and will eventually be completely replaced with a more general post-processor. 

Several new problem specific routines were added. A routine to update boundary values during a run was added. 
A number of other changes and bug-fixes were also made, primarily involving geometry specific bugs. Also, the build 
process was updated, relying on improved makefiles and a custom Perl script. 

To insure that ZEUS-MP is portable, we have tested and run benchmarks on several different systems. We are 
aided by the fact that ZEUS uses NCSA's hdf4 format as its primary output format. By using a portable output 
format, we do not have to worry about endian issues with our output files. The one place we break this is in the 
restart dumps which are produced as unformatted FORTRAN binary data. We have run and used ZEUS-MP on a 
variety of GNU/Linux distributions (both 2.4.x and 2.6.x kernels) on both workstations and a Beowulf cluster. AMD 
Athlon processors, AMD Opteron processors (in 32 bit mode only as that is all we currently have available), and Intel 
Pentium 4 processors have primarily been used. We have also compiled and run ZEUS-MP on Apple's Mac OS X 
(darwin) on the G4 PowerPC processor. We have also attempted to run ZEUS-MP on a Sun Ultra80 (UltraSPARC 
II processor) running Solaris 8.0. As we were not able to successfully use the MPI library needed, we were not able 
to run ZEUS-MP. This is the one weakness in our portability. If a working MPI library (preferably lam-mpi) and the 
hdf4 library cannot be compiled, ZEUS-MP cannot run on a system. This should not be a problem on any modern, 
Unix-like system. 

A.l. Performance 

To test performance and scaling, a simulation was run with a 100 3 Cartesian grid for 200 time steps in pure 
hydrodynamics mode with radiative cooling added. This simulation was run as one to eight processes (not necessarily 
processors). For single processor machines, this means multiple processes on one processor (and presumably a large 
performance penalty). For the cluster tests this is not a problem. The timing results are shown in Figure lAT2l There 
are two surprises in these results. The first is that at 8 processors, code compiled with the GNU compilers seems to 
outperform the code from the Intel compiler (when used on AMD processors). This is surprising as the Intel compiler 
is widely considered to produce the fastest code for x86 processors. The second surprising feature is that the G4, while 
it performs far below anything else in the tests, does not appear to have any penalty for multiple processes on a single 
processor. 

The results from the scaling test are shown in Figure IA13I While at no point beyond two processors do we get the 
theoretical linear increase, we continue to get a decent speedup (with a constant slope in most cases) up to at least 8 
processors, with the code produced with the GNU compiler getting much closer to linear than the Intel compiler. 
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